Quartically regularized polynomials

Adapted from: Section 3.3.1 of (Cartis and Zhu, 2026)

Higher-order optimization methods minimize, at each iteration, a cubic Taylor model regularized by a quartic term:

\[m(s) = f_0 + g^\top s + \frac{1}{2} H[s]^2 + \frac{1}{6} T[s]^3 + \frac{\sigma}{4} \|s\|_2^4.\]

For any $\sigma > 0$, this model is bounded below so $m(s) - m^*$ is nonnegative where $m^*$ is its global minimum. If $m(s) - m^*$ is moreover a sum of squares then $m^*$ can be computed exactly with a single Sum-of-Squares program, i.e., the zeroth level of the Lasserre hierarchy. (Cartis and Zhu, 2026) shows that this is the case for $\sigma$ large enough and that, when $T = 0$, it is the case for all $\sigma > 0$. They also show that the geometry of the regularization matters: with the separable quartic regularization $\|s\|_4^4 = s_1^4 + \cdots + s_n^4$ instead of the Euclidean one $\|s\|_2^4 = (s_1^2 + \cdots + s_n^2)^2$, $m(s) - m^*$ may fail to be a sum of squares for every $\sigma > 0$. In this tutorial, we verify this numerically on their example.

using LinearAlgebra
using DynamicPolynomials
@polyvar s[1:3]

(DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}[s₁, s₂, s₃],)

We consider the following nonconvex quadratic part whose Hessian has eigenvalues $-4$, $-4$ and $20$:

quad = 2 * sum(s .^ 2) + 8 * (s[1] * s[2] + s[1] * s[3] + s[2] * s[3])

\[ 2s_{3}^{2} + 8s_{2}s_{3} + 2s_{2}^{2} + 8s_{1}s_{3} + 8s_{1}s_{2} + 2s_{1}^{2} \]

We need to pick an SDP solver, see here for a list of the available choices.

using SumOfSquares
import Clarabel
solver = optimizer_with_attributes(Clarabel.Optimizer, MOI.Silent() => true)

MathOptInterface.OptimizerWithAttributes(Clarabel.MOIwrapper.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.Silent() => true])

The following function computes the largest $\gamma$ such that multiplier * (p - γ) is a sum of squares. With the default multiplier, this is the zeroth level of the Lasserre hierarchy. It returns the lower bound $\gamma$ and the constraint reference, which we will use to extract minimizers from the moment matrix.

function sos_lower_bound(p; multiplier = 1)
    model = SOSModel(solver)
    @variable(model, γ)
    @objective(model, Max, γ)
    con_ref = @constraint(model, multiplier * (p - γ) >= 0)
    optimize!(model)
    @assert termination_status(model) == MOI.OPTIMAL
    return value(γ), con_ref
end

sos_lower_bound (generic function with 1 method)

Separable quartic regularization

We first consider the separable regularization with $\sigma = 4$. The resulting polynomial was shown in (Ahmadi et al., 2023; Theorem 3.3) to be such that $m_A(s) - m_A^*$ is not a sum of squares.

m_A = quad + sum(s .^ 4)

\[ 2s_{3}^{2} + 8s_{2}s_{3} + 2s_{2}^{2} + 8s_{1}s_{3} + 8s_{1}s_{2} + 2s_{1}^{2} + s_{3}^{4} + s_{2}^{4} + s_{1}^{4} \]

The zeroth level of the Lasserre hierarchy gives the following lower bound:

γ_A, _ = sos_lower_bound(m_A)
γ_A

-2.385952772459071

Is this lower bound tight ? Recall that a polynomial $p$ is nonnegative if $(s_1^2 + s_2^2 + s_3^2) p(s)$ is a sum of squares, see the Motzkin tutorial. Using this multiplier, we get a larger lower bound:

γ_A_mult, con_ref = sos_lower_bound(m_A, multiplier = sum(s .^ 2))
γ_A_mult

-2.1442784229172616

We can extract the minimizers from the moment matrix:

atoms = atomic_measure(moment_matrix(con_ref), 1e-4)

Atomic measure on the variables s[1], s[2], s[3] with 7 atoms:
 at [0.6674021396193536, -1.1497297911622102, 0.6673976689402396] with weight 0.07442764536865731
 at [-0.6674021522297394, -0.6673980177637419, 1.1497302302674677] with weight 0.07044123967344412
 at [1.1497283431423126, -0.667399018846945, -0.6673987145205218] with weight 0.08109682421963398
 at [-1.1497283004340575, 0.6673990116055337, 0.6673986677456719] with weight 0.08109683782196779
 at [-0.6674012443392733, 1.1497294672965874, -0.6673979692274824] with weight 0.07442773228171859
 at [0.6674012142414039, 0.6673983488707749, -1.1497298832053748] with weight 0.07044132634156423
 at [0.0, 0.0, 0.0] with weight 0.0

We recover the six global minimizers $\pm(-1.1498, 0.6674, 0.6674)$ and their permutations. The moment matrix also contains an atom at the origin with zero weight. It comes from the multiplier which vanishes at the origin, so we ignore it.

minimizers = [atom.center for atom in atoms.atoms if atom.weight > 1e-4]

6-element Vector{Vector{Float64}}:
 [0.6674021396193536, -1.1497297911622102, 0.6673976689402396]
 [-0.6674021522297394, -0.6673980177637419, 1.1497302302674677]
 [1.1497283431423126, -0.667399018846945, -0.6673987145205218]
 [-1.1497283004340575, 0.6673990116055337, 0.6673986677456719]
 [-0.6674012443392733, 1.1497294672965874, -0.6673979692274824]
 [0.6674012142414039, 0.6673983488707749, -1.1497298832053748]

Evaluating $m_A$ at the minimizers confirms that this is the global minimum and that the zeroth level of the hierarchy has a gap:

m_A_min = minimum(x -> m_A(s => x), minimizers)
m_A_min

-2.1442784126213783

As shown in (Cartis and Zhu, 2026; Lemma 3.3), the gap does not close as $\sigma$ increases. Indeed, with $\tilde{s} = \sqrt{\sigma}s/2$, we have $\mathrm{quad}(s) + \sigma \|s\|_4^4/4 = 4 m_A(\tilde{s})/\sigma$ so the ratio between the lower bound and the minimum is independent of $\sigma$:

for σ in [1, 10, 100]
    m = quad + σ / 4 * sum(s .^ 4)
    γ, _ = sos_lower_bound(m)
    println("σ = $σ: lower bound = $γ, minimum = $(4m_A_min / σ)")
end

σ = 1: lower bound = -9.543811171212013, minimum = -8.577113650485513
σ = 10: lower bound = -0.95438112377011, minimum = -0.8577113650485513
σ = 100: lower bound = -0.09543811863380461, minimum = -0.08577113650485513

Euclidean quartic regularization

Let us now replace the separable regularization with the Euclidean one:

σ = 4
m_E = quad + σ / 4 * sum(s .^ 2)^2

\[ 2.0s_{3}^{2} + 8.0s_{2}s_{3} + 2.0s_{2}^{2} + 8.0s_{1}s_{3} + 8.0s_{1}s_{2} + 2.0s_{1}^{2} + s_{3}^{4} + 2.0s_{2}^{2}s_{3}^{2} + s_{2}^{4} + 2.0s_{1}^{2}s_{3}^{2} + 2.0s_{1}^{2}s_{2}^{2} + s_{1}^{4} \]

This is a quadratic model with Euclidean quartic regularization ($T = 0$) so $m_E(s) - m_E^*$ is a sum of squares for every $\sigma > 0$. Indeed, it can be checked that

\[m_E(s) + \frac{4}{\sigma} = \frac{\sigma}{4}\left(\|s\|_2^2 - \frac{4}{\sigma}\right)^2 + 4(s_1 + s_2 + s_3)^2.\]

The zeroth level of the hierarchy finds this lower bound $-4/\sigma = -1$:

γ_E, con_ref = sos_lower_bound(m_E)
γ_E

-0.9999999989936235

The global minimizers form the circle $s_1 + s_2 + s_3 = 0$, $\|s\|_2 = 2/\sqrt{\sigma}$. Even if there are infinitely many of them, the low rank decomposition of the moment matrix still finds a few points of this circle:

atoms = atomic_measure(moment_matrix(con_ref), 1e-4)

Atomic measure on the variables s[1], s[2], s[3] with 4 atoms:
 at [0.5773553649081746, 0.2113282386431019, -0.7886836035398079] with weight 0.24601237796659883
 at [-0.577355365380753, 0.7886836038617806, -0.21132823855682908] with weight 0.2460143763033419
 at [-0.5773553649081774, -0.21132823864310085, 0.7886836035398095] with weight 0.24601237796659872
 at [0.5773553653807553, -0.7886836038617819, 0.211328238556828] with weight 0.2460143763033413

We can verify that they belong to the circle:

for atom in atoms.atoms
    x = atom.center
    println("sum = $(sum(x)), norm = $(norm(x))")
end

sum = 1.1468603844377867e-11, norm = 1.0000103341109952
sum = -7.580153171815596e-11, norm = 1.0000103346195384
sum = -1.1468825888982792e-11, norm = 1.0000103341109978
sum = 7.580133742912665e-11, norm = 1.0000103346195404

Adding a cubic term

We now add a small cubic term $0.1(s_1^3 + s_2^3 + s_3^3)$ to both models. With the separable regularization, there is still a gap between the lower bound of the zeroth level and the global minimum:

cubic = 0.1 * sum(s .^ 3)
m_A2 = m_A + cubic
γ_A2, _ = sos_lower_bound(m_A2)
γ_A2

-2.3989905400563516

The multiplier closes the gap again and gives the global minimum:

γ_A2_mult, con_ref = sos_lower_bound(m_A2, multiplier = sum(s .^ 2))
γ_A2_mult

-2.240948485668829

With the Euclidean regularization, the model is now a cubic model so (Cartis and Zhu, 2026; Theorem 2.3) only guarantees that $m(s) - m^*$ is a sum of squares for $\sigma$ large enough. We compare the lower bound found by the zeroth level of the hierarchy with the value of $m$ at the minimizers extracted from the moment matrix. If they match, then the lower bound is the global minimum.

for σ in [0.1, 1, 4, 10, 40]
    m = quad + cubic + σ / 4 * sum(s .^ 2)^2
    γ, c = sos_lower_bound(m)
    μ = atomic_measure(moment_matrix(c), 1e-4)
    m_min = minimum(atom -> m(s => atom.center), μ.atoms)
    println("σ = $σ: lower bound = $γ, value at extracted minimizers = $m_min")
end

σ = 0.1: lower bound = -53.43765171837269, value at extracted minimizers = -53.437653996586846
σ = 1.0: lower bound = -4.353277438109152, value at extracted minimizers = -4.353277564505868
σ = 4.0: lower bound = -1.042438226708667, value at extracted minimizers = -1.0424382384183706
σ = 10.0: lower bound = -0.4105830614272902, value at extracted minimizers = -0.4105830511102412
σ = 40.0: lower bound = -0.10130681226787011, value at extracted minimizers = -0.10130677616184373

For this example, the zeroth level is exact even for small values of $\sigma$. The theoretical threshold of (Cartis and Zhu, 2026; Theorem 2.3) is sufficient but not necessary.

This page was generated using Literate.jl.